Cuticular hydrocarbon profiles reveal geographic chemotypes in stingless bees (Hymenoptera: Meliponini)

Cuticular hydrocarbon (CHCs) variation has been detected in various insect taxa, but the potential contribution of cuticular compounds for analyzing intraspecific diversity at the population level has been little explored. Here we assess for the first time intraspecific variation in the CHC profile of stingless bees, using the species Melipona beecheii and Nannotrigona perilampoides. The objective is determining whether intraspecific variation can be useful for population identification. We found species-specific chemical patterns and extensive variation within each species. Notably, chemotypes were significantly associated to geographic origin in N. perilampoides but less so in M. beecheii and we discuss possible explanations for these patterns. Our results support the use of CHCs in conjunction with other methods in emerging problems such as undetected colony mobilization across regions. As CHCs are involved in several aspects of stingless bee recognition and interactions, it would be essential to unravel how these chemical signatures evolve across populations.

The composition and quantity of compounds on the surface of the cuticle of insects can be species-specific 1 .Thus, differences in cuticular profiles have been used as valid tools for species identification across several insect groups and for closely related ones 2 .As taxonomic tools, cuticular profiles are successfully combined with morphometrics and molecular tools for delimiting species 3,4 , including cryptic ones 5,6 .Notably, although chemotypes have been detected in various insect taxa [6][7][8] , the potential contribution of cuticular compounds for analyzing intraspecific diversity has been little explored.
The majority of insect cuticular compounds are hydrocarbons (CHCs), mainly alkanes, alkenes, and branched alkanes, with a chain length typically varying from C17 to C35, that are synthesized in specialized cells called oenocytes 9,10 .CHCs are subsequently transported (together with fatty acids, alcohols, glycerides, phospholipids and glycolipids) to the cuticle of the insect where they form a lipid wax layer that plays critical roles in the integrity and communication of individuals 1,11 .Primarily, CHCs serve as barrier to avoid water loss and abrasion of the insect's body and protection against microorganisms 3,10 .In addition, in social insects CHCs can also convey information about the reproductive and physiological status, hierarchy, and task activity of different individuals [12][13][14][15][16][17] .
Most intraspecific-studies of CHCs in eusocial bees, including the Meliponini (Anthophila: Apidae) also known as stingless bees, have centered around differences among individuals and their role in communication.Thus, variation in cuticular profiles has been observed when comparing castes (queens and workers) in Melipona bicolor 18 or among workers of different age and tasks groups in Frieseomelitta varia 19 , Melipona marginata 20 , Friesella schrottkyi 21 , Tetragonisca angustula 22 and Melipona solani 23 .Although the importance of these compounds in the organization and function of colonies, is clear, considerably less is known about within species variation at the population level 3,24 .
Mesoamerica is a region where several stingless bee species have been cultured over centuries.Notably, the number of species of these insects under husbandry is increasing locally and worldwide 25 .Due to their popularity and escalating demand, translocation of stingless bee colonies across geographic regions is more frequently occurring 26,27 .Translocating colonies across regions or to regions naturally lacking stingless bees,

Chemical analysis of extracts
Extracts were analyzed on an Agilent Technologies 7890A gas chromatograph coupled to a 5975C Mass spectrometer.The separation was carried out with an Agilent HP-5MS capillary column (30 m × 250 µm × 0.25 µm).The oven was programmed with a temperature ramp from 40 to 240 °C at a speed of 40 °C/min, held for 5 min and then increased to 300 °C at 10 °C/min and held for 9 min.Helium was used as a carrier gas with a constant flow rate of 1 mL/min.For each analysis, 2 µL of wings extract was injected into the GC inlet at 300 °C in splitless mode (100 mL/min at 1.5 min).The mass spectrometer operated in electron impact mode at 70 eV.The temperature of the ion source was 230 °C and 150 °C for the quadrupole.A series of chromatograms were studied with the Agilent ChemStation software 33 in which the compounds of the cuticular extracts were identified and characterized.

Characterization and quantification of cuticular compounds
A standard mixture of C7 to C40 hydrocarbons (Merck, 49,452-U) with 5 µg/mL of each compound in hexane was analyzed to obtain their retention times and calculate the Kovats retention index of the unknown compounds.The retention times and indices obtained were subsequently compared with those reported in the literature or standard compounds for compound identification (see Supplementary files Table S1 online).The mass spectral patterns were also analyzed and compared to the NIST MS library 34 .The experimental Kovats retention indices (KI) were calculated with the formula: The n symbolizes the number of carbons of the n-alkane that eluted before the unknown compound, t x the retention time of the unknown compound, t n the retention time of the n-alkane that eluted before t x and t n+1 the time of the one that eluted immediately after t x .
The relative proportion of each identified CHC was calculated using the quotient of the sum of the peak areas of a given compound divided by the overall sum of the peak areas of all compounds for each locality.From the overall proportion of each compound, we calculated the proportion of alkanes and alkenes for each species.A general comparison of chemical profiles between M. beecheii and N. perilampoides was done by plotting the proportion of the main compound categories of CHCs, namely, alkanes, alkenes and their isomers.The relative proportion percentage (%) of each CHC in relation to the total composition of the extracts from each location was calculated.Comparison among geographic regions for each CHC per species were done by means of the Kruskal-Wallis statistic with the function 'kruskal.test' 35followed by Dunn's multiple comparison tests and the corrected Bonferroni values, using using the ' dunnTest' function of the 'FSA' package 36 .

Comparison of chemical profiles
To compare the chemical profiles among geographic regions within each species, a non-metric multidimensional scaling (NMDS) analysis was performed using the 'nmds' function in R package 'ecodist' 37,38 .The NMDS analysis was conducted using a triangular matrix generated as of the relative abundances.Based on the triangular matrix, we built a two-dimensional NMDS plot using 50 iterations per run in the R package 'ecodist' 38 .To ensure the reproducibility of the NMDS plot, we used the 'set.seed'function 35 .In the NMDS plot, deviations are expressed in terms of "stress" where values below 0.15 indicate a good fit for the data 5 .Differences between regions were evaluated by one-way analysis of similarity (ANOSIM); both NMDS and ANOSIM were based on the Bray-Curtis dissimilarity index implemented in the R package 'vegan' 39 .ANOSIM test is a special form of

Results
The cuticular hydrocarbon profile of both studied species was contrasting.Despite sharing various hydrocarbons, the relative frequency of each compound significantly varied between species.It is also evident that the profile of M. beecheii included isomers of some alkenes, but in contrast N. perilampoides lacked isomers in its profile  www.nature.com/scientificreports/(Fig. 2).It is therefore, as expected, possible to differentiate both species in accordance to their general profile of cuticular hydrocarbons.The chemical analysis of the cuticle of M. beecheii identified a total of 20 CHCs, three of them with isomers (Fig. 2A, Table 4).In the case of N. perilampoides the analysis of the cuticular profile identified a total of 16 CHCs, none of them with isomers (Fig. 2B, Table 5).
In Tables 2 and 3 we present the average relative proportion for the compounds found in the typical cuticular profile for colonies of M. beecheii and N. perilampoides, respectively, in each geographical location sampled.
For M. beecheii the chain length of the identified CHCs varied between C19 and C33, which were classified as linear alkanes and alkenes.The most frequent CHCs in the populations studied were C25, C27, C29, C31 and C33 with their respective alkenes.Notably, the proportion of each of these CHCs varied across populations (Fig. 3A).The populations from PY and Tab East had less proportion of the alkane C25 and alkenes C25:1-1 and C25:1-2, as well as the alkane C27 and its alkenes C27:1-1 and C27:1-2, both of which were more frequent in the other four populations.In contrast, PY and Tab East had significantly higher proportion of C31 and its alkene C31:1 compared to the other regions.Likewise, alkane C33 was mostly detected in PY and Tab East compared with the other four regions where its presence was negligible (Fig. 3A).Thus, alkanes and alkenes of larger chain molecules were found in PY and Tab East, while CHCs with shorter chain lengths were found in the other four regions.
In the case of N. perilampoides, the CHCs profile of this species varied between C21 and C31, which were also classified as linear alkanes and alkenes (Table 3).Overall, in N. perilampoides, it was more evident that the profile included a larger proportion of alkenes compared with alkanes, compared with M. beecheii (Fig. 4A).The most abundant CHC was the alkene C29:1 but its alkane was found in small proportion and in Tab was practically not detected.Other most frequent CHCs in the populations studied were C25 and C27 with their respective alkenes.
The plot of NMDS for M. beecheii showed clear separation among the different populations in accordance to their geographic origin (Fig. 3B).The results of the ANOSIM further supported this finding (Table 3).The overall value of the ANOSIM for M. beecheii was 0.859.The most distant populations regarding their hydrocarbon profiles were Tab East and Chia (ANOSIM = 1), followed by Tab West and Chia (ANOSIM = 0.978) and Tab East with Pue (ANOSIM = 0.975).Other ANOSIM values showing strong differences among paired populations (> 0.9) were those of Pue with Chia, PY with Tab West and the latter with Tab East.No significant differences were found between the paired populations of Tab West with Ver, Ver and Oax and Oax with Pue (this one had the lowest ANOSIM value = 0.038).The results of the UPGMA confirmed that it is possible to identify three main chemotypes among the studied populations of M. beecheii (Fig. 4A).One group formed by the bees from PY and East Tab, a second larger group formed by Tab West, Ver, Oax and Pue and finally a third group constituted solely by the bees of Chia alone.Some aspects worth to note are first, that the populations of M. beecheii from the state of Tabasco (Tab) can be clearly separated into two chemotypes, one from the East whose profiles are similar to the bees from PY and another of the Tab West bees which are more similar in their cuticular profiles to the bees from Ver, Pue and Oax.Second, the population of Chia forms a very compact and clearly separated group from all others.Lastly, two colonies from PY were placed close to the population from Oax, this is interesting, and it would be important to verify the origin and the possibility of colony exchange among these regions.
Similar to M. beecheii, the proportion of the different CHCs varied across populations of N. perilampoides (Fig. 5A).The populations from PY and Tab had the highest proportion of the alkene C29:1.Interestingly, this alkene seems to have a decreasing proportional gradient from the populations of East of Mexico towards the West, with the population on the Pacific state of Jal having the lowest proportion of this compound.A reverse situation seems to occur with alkane C25 which contributed proportionally less to the profile of populations in PY and Tab compared with the other four (Fig. 5A).Notably, the more frequent compound found in the profile of N. perilampoides from Jal was the alkene C25:1.Thus, alkenes of larger chain molecules were also found in N. perilampoides from PY and Tab, while compounds with shorter chain molecules were found in the other four regions.
The NMDS plot for N. perilampoides showed separation among the different populations in accordance to their geographic origin (Fig. 5B).However, the value of the overall ANOSIM for this species was comparatively lower = 0.603 (Table 4), meaning that the separation among some populations was less pronounced.The www.nature.com/scientificreports/most distant populations of N. perilampoides were Pue and Jal (ANOSIM = 0.996), followed by Jal and Oax (ANOSIM = 0.96) and Jal with PY (ANOSIM = 0.933).Other ANOSIM value showing strong differences among paired populations (= 0.868) was that of Pue with PY.Notably, six comparisons among paired populations showed no significant differences for N. perilampoides.No significant differences were found between Tab with Ver, Tab and Oax, Tab and Pue, Ver and Oax, Oax with Pue and Ver with Pue (this paired comparison had the lowest ANOSIM value = − 0.102).Given these results, it is also possible to identify three main chemotypes in the UPGMA of the studied populations of N. perilampoides (Fig. 6).Clearly different to all regions is the population from Jal on the Pacific Coast, as well as the PY population which was significantly different to all others.A third larger group is formed by Tab, Ver, Oax and Pue (Fig. 6A).One important aspect to note is the negative value obtained for the paired comparison of Ver and Pue (Table 4 and Fig. 6A).This pattern usually results when the difference within one population is larger than the difference between populations 42 .In fact, on the NDSM plot it is possible noticing that while some colonies from Ver are similar to Pue and Oax, some colonies from Ver cluster separately from the former (Fig. 5B).www.nature.com/scientificreports/ The results of the Mantel test correlating chemical with geographic distance were highly significant in N. perilampoides (r M = 0.804, P = 0.006) , as can be also seen on Fig. 6A and B. This result indicates isolation by distance as a possible factor causing differentiation of the populations in N. perilampoides.However, the value of the Mantel test was not significant in M. beecheii (r M = 0.33, P = 0.094) , indicating that geographic separation is a comparatively less powerful factor in the differentiation among populations of this species (Table 5).

Discussion
In general, there is still a lack of solid evidence on the natural intra-species variation in CHC profiles, even in well-established insect species 3 and these compounds have thus, not been extensively used for chemotaxonomy in social hymenopterans 43 .We analyzed the profile of CHCs of the Mesoamerican stingless bees M. beecheii and N. perilampoides and found clear species-specific chemical patterns.Our results revealed extensive variation within each species and notably, that chemotypes were significantly linked to the geographic regions from where populations originated.This evidence supports the use of CHCs as additional markers for the study of variation and evolution in stingless bees.
Normally, CHCs vary between species in the types of compounds that they produce and/or the relative quantities of shared compounds 3,4 .Intraspecific CHC variation in the closely related ants and wasps is mostly quantitative, that is, individuals possess the same set of hydrocarbons, but in different relative quantities 44,45 .Likewise, our study in stingless bees revealed that in the two species studied, variation in the cuticular profiles could be found between species but also within species.Variation was mostly quantitative, but the profiles obtained even for similar CHCs was sufficient to clearly separate both species and to identify geographic populations within them.
Noteworthy, chemotypes were detected within each species and these were significantly determined by the geographic region.We found that alkanes with shorter chain lengths (molecules with a comparatively higher fluidity) were more frequent in N. perilampoides compared with M. beecheii.Moreover, more compounds were detected for both species in the regions with more constant humidity all year round compared with the Yucatan Peninsula which experiences a period of drought between March and May.Indeed, the Yucatán Peninsula is a separate biogeographic region with a distinct climate, predominantly of the Aw type with comparatively less rainfall (subhumid).In contrast, Veracruz and Tabasco are predominantly included in the Af and Am types, characterized by average higher rainfall 46 .Thus, the particular profile seems to be in concordance to the habitats in which both species are naturally present.Warm and dry conditions can lead to higher proportions of more aggregating substances like n-alkanes, while the more fluid unsaturated hydrocarbons tend to decrease 45 .N. perilampoides is a bee that can be found from altitudes at sea level up to 1500 m, even reaching cloud tropical forests where temperatures can be relatively low while M. beecheii is usually found at sea level not exceeding a couple hundred meters of altitude 30 .Even more, the dry conditions of the lowlands of Yucatan seem to also affect the pattern of CHCs with less diverse compounds found compared with the other regions.Such differences in habitat preference by both species and of the geographic region where their populations are found could be determinant in the pattern of cuticular hydrocarbons found in each.Stress value for NMDS configuration = 0.07.Acronyms are those in Table 1.www.nature.com/scientificreports/Compounds within each species were identified that could be useful in the identification of geographic populations.In the case of M. beecheii, the alkene C31:1 was present in frequencies of 35% and above in the Yucatan Peninsula and Tab East.This compound was only marginal present in all the other populations, except in Oaxaca where it was present but with a frequency below 15%.Likewise, the frequency of C25 and C27 and its respective alkenes, were only marginal in populations from the Yucatan Peninsula and Tabasco East, but it was more readily identified in the other four populations, with the exception of Chiapas where C25 was also present at low frequency.Noteworthy, the frequency of C27 and C27:1 was highest in Chiapas where both represented more than 50% of the cuticular profile.In the case of N. perilampoides the alkene C29:1 was present at frequencies of 35% and above in populations from the Yucatan Peninsula and Tabasco, but below these levels in all other populations.The alkane C25 was also more frequent in populations other than the Yucatan Peninsula and Tabasco where its frequency was below 20%.Interestingly, the population from Jalisco can be readily identified by the high levels of the alkene C25:1 where it was the most abundant compound (above 40%) compared with all other populations where it was present at very low frequencies.
Geographical barriers, coupled with related ecological differences, may have driven population divergence and possibly explain the emergence of chemical differences and geographic chemotypes.In addition, orographic and biotopic particularities have probably acted selecting for particular CHCs to protect individuals under different environmental conditions, thus, favoring differentiation even further 47 .Generally, stingless bees show strong genetic structure [48][49][50][51] .This is a consequence of their particular mode of reproduction related to dependent colony foundation, in which a virgin queen swarms with a group of workers to start a new colony 52 .However, in contrast with other highly eusocial insects, newly formed colonies of stingless bees maintain long-time bonds for food and resources with the mother colony 25,53 .The link with the mother colony reduces the spreading of new colonies (and reproductive females) over long spaces.This favors genetic differentiation that is frequently associated with geographic distance in various species of stingless bees 47,48,54 .
Additionally, differences in chemical profiles may be one feature of the adaptation of populations to specific ecosystems, an essential aspect when considering colony trade.This implies that some species might suffer and perish due to intolerance to new climatic features 26 .In this regard, one practical aspect of our work is the identification of populations of two stingless bee species whose colonies are increasingly commercialized with a  1 and data are from Table 3. www.nature.com/scientificreports/latent risk of being shipped to distant areas, many outside their original geographic region 27 .The effect of man induced translocation is evident by the weaker genetic structure of managed species of stingless bees compared to feral ones 51,55 .In Brazil for instance, the translocation of colonies poses a severe problem for species genetic integrity.At least 33 species of stingless bee have been negotiated and over 40% of sellers and buyers were located outside the natural range of the species 26,56 .Such trade has mixed populations and introduced species, allows the potential dissemination of hitchhiker symbionts and diseases, and mostly does not follow national legislation policies 26,56 .Interestingly, the results from the Mantel test indicate that N. perilampoides populations can be more clearly separated in accordance with the geographic region compared with M. beecheii.This could be an indication that populations of N. perilampides seem more structured, and probably no admixture has occurred as consequence of colony translocation for this species.In contrast, differences were less pronounced in accordance to geographic distance in M. beecheii.It may be possible that given the significance of this species for ancient stingless beekeeping, translocation of some colonies may have occurred in the past 25 .However, this deserves further analysis, possibly including other markers and more colonies, as it may also be an indication of recent admixture among populations.With the increasing popularity of stingless bees, in particular of M. beecheii, it would be important to monitor the movement of colonies across regions and CHCs could provide a rapid and reliable tool in conjunction with other methods 4,48,57 .
CHCs have high plasticity, allowing insects to quickly adapt their chemical profiles to external selection pressures 58 .Because of this, their use as taxonomic tools can be important in solving problems where other methods are not capable of detecting differences.In a pair of sibling species of orchid bees, DNA barcoding based on changes in the mitochondrial cytochrome c oxidase subunit 1 (COX1) showed no significant differentiation 59 , while microsatellite markers showed too much overlap in allele sizes to be diagnostic 60 .The species could nevertheless be distinguished based on their CHCs profiles 5 .Recently, geographic populations of N. perilampoides, showed no differentiation after the analyses of the nuclear regions ITS1, ITS2 nor the mitochondrial COX1.Differences were only observed for the mitochondrial 16S marker.However, different haplotypes differed only in one or two mutational steps, so in general, the nucleotide diversity was low for this marker 48 , rendering molecular tools of scarce utility for analyzing differentiation among populations of this species.However, as our study showed, significant differences could be found in the profile of CHCs of geographic populations of this species, rendering CHCs as a good set of epigenetic markers.
One important aspect for the use of CHCs as taxonomic tools is the genetic basis behind their variation among species 61 .In ants, recent studies have shown that cuticular profiles are heritable and species-specific 62,63 .In stingless bees, alkanes, alkenes and alkadienes showed no or little correlation with the geographical distribution of different species 64 .Similarly, chemical differentiation of thirteen species of Melipona bees showed no significant effect of environmental conditions 32 .These results indicate strong stability of CHC profiles over large geographical distances among species.Nevertheless, practically nothing is known about the relative contribution of environmental and genetic factors in determining intraspecific differentiation in stingless bees (and other social insects) making this is a pending and essential aspect to further study.
Our study found intraspecific CHCs variation in stingless bees which support their use as chemotaxonomic tools for the analysis of geographical differences within species.CHCs data bases could be included in automated systems in conjunction with other methods for an accurate identification of stingless bee populations and trace possible movements in the future.The relative influence of genetic and environmentally factors explaining CHC profile variation are still not clear and deserve further analysis.As CHCs are involved in several aspects of bee recognition and interactions, it would be key to unravel how these chemical signatures evolve across populations 45,65 .

Figure 1 .
Figure 1.Map showing the localities where samples of M. beecheii (A) and N. perilampoides (B) were collected in different geographic regions across Mexico.Names of localities and acronyms for regions are those in Table1.

Figure 3 .
Figure 3. Relative frequencies of alkanes (solid bars) and alkenes (pattern bars) in the cuticular profile of M. beecheii colonies collected in different states and geographic regions across Mexico (A).In (B), Nonmetric multidimensional scaling (NMDS) plot of cuticular hydrocarbon variation among populations of M. beecheii.Stress value for NMDS configuration = 0.07.Acronyms are those in Table1.

Figure 4 .
Figure 4. UPGMAs for M. beecheii, in (A) based on chemical distances and (B) based on geographic distance among different regions of Mexico.Acronyms are those in Table1and data are from Table3.

Figure 6 .
Figure 6.UPGMAs for N. perilampoides, in (A) based on chemical distances and (B) based on geographic distance among different regions of Mexico.Acronyms are those in Table1and data are from Table3.

Table 1 .
Number of colonies of M. beecheii and N. perilampoides sampled per Region, State, Locality and meliponary.*All colonies from this locality were feral.Capital letters in brackets within Region denote the acronym used to refer to them in the text and in subsequent Tables and Figures.

Table 2 .
Relative proportion (± Standard error) of cuticular hydrocarbons found in M. beecheii for different states and regions of Mexico.The value of the Kruskal-Wallis test is represented by the letter H.The values in parenthesis represent the average ranks derived from Dunn's multiple comparison test.Different letters denote significant differences at p < 0.05.

Table 3 .
Relative proportion (± Standard error) of cuticular hydrocarbons found in N. perilampoides for different states and regions of Mexico.The value of the Kruskal-Wallis test is represented by the letter H.The values in parenthesis represent the average ranks derived from Dunn's multiple comparison test.Different letters denote significant differences at p < 0.05.

Table 5 .
Values for paired comparison of chemical distance of N. perilampoides colonies from different geographic regions (values closer to 1 indicate larger differentiation).The global value of ANOSIM for all groups was R = 0.603.Figures in bold and italic indicate non-significant differences between populations (p > 0.01).Below the diagonal geographic distance between regions given in km.Acronyms for regions are those in Table1.